############################## bkmrbkmr总体效应 #############################
##### 科研工具7   bkmrbkmr总体效应
library(ggplot2)
library(bkmr)
library(ggplot2)
library(readr)
setwd('/home/r_temp//202443/1/7/')

demo <- read_csv("demo.csv")
# demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
#                  'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
#                  'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
#                  'Hyperuricemiaint','LBDSUASI')]

tryCatch({
demoS4 <- demo[,c('MCOP','MnBP','MiNP','MECPP','Hyperuricemiaint')]
},error = function(e){
print("获取数据出错")
})

if(dim(demoS4)[1]>=300){
demoS4 <- demoS4[c(1:300),]
}

# demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
# Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
#                'MCOP','MCPP','MEP','MBzP')]

tryCatch({
Z <- demoS4[,c('MCOP','MnBP','MiNP','MECPP')]
},error = function(e){
print("获取数据出错")
})



tryCatch({
for(i in 1:4){Z[,i] <- log(Z[,i])}
},error = function(e){
print("log 转换出错提示")
})
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 200, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
risks.overall = OverallRiskSummaries(fit=fitkm,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="OverallRiskSummaries.tiff", width=8, height=8,units="in",dpi=150)
ggplot(risks.overall,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
geom_hline(yintercept = 0,lty=2,col='#BA0707')+
geom_pointrange()
dev.off()